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Abstract 

A cell is polarised when it has developed a main axis of organisa- 
tion through the reorganisation of its cytosqueleton and its intracellular 
organelles. Polarisation can occur spontaneously or be triggered by ex- 
ternal signals, like gradients of signaling molecules ... Following [5] and 
[3], in this work, we study mathematical models for cell polarisation. 
These models are based on nonlinear convection-diffusion equations. 
The nonlinearity in the transport term expresses the positive loop be- 
tween the level of protein concentration localised in a small area of the 
cell membrane and the number of new proteins that will be convected 
to the same area. We perform numerical simulations and we illustrate 
that these models are rich enough to describe the apparition of a po- 
larisome. 

Keywords: Cell polarisation, global existence, blow-up, numerical 
simulations, Keller-Segel system. 

1 Introduction 

Cell polarisation is a major step involved in several important cellular pro- 
cesses such as directional migration, growth, oriented secretion, cell division, 
mating or morphogenesis. When a cell is not polarised molecular markers 
(proteins CDC42) are uniformly distributed on the membrane while polari- 
sation is characterized by the concentration of molecular markers in a small 
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area of the cell membrane. In [T2], it has been observed that if the external 
pheromone concentration is above a critical concentration, polarisation can 
occur spontaneously. It has also been observed that cell asymmetry can be 
driven by an external asymmetric stimulus. 

Cell polarisation in yeast cells has been intensively studied during the 
past decade. Recently, many models describing cell polarisation have been 
developed. The majority of these models are based on reaction-diffusion 
systems where polarisation appears as a type of Turing instability [9] , [12] , 
[TO] , or due to stochastic fluctuations [2], other models include cytoskele- 
ton proteins as a regulatory factor [6], [H]. Many biological studies have 
shown that the cytoskeleton plays an important role in polarisation. It has 
been suggested that the cytoskeleton has a positive feedback on molecu- 
lar markers density. Indeed, disruption of transport along the cytoskeleton 
greatly reduces the stability of polar cap [H]. The cell cytoskeleton is a 
network of long semi-flexible filaments made up of protein subunits [13] . 
These filaments (mainly actin or microtubules) act as roads along which 
motor proteins are able to perform a biased ballistic motion and carry vari- 
ous molecules. Molecular markers play a key role in the formation of these 
filaments. 

Following [8], [1] and [3], in this work we study models that describe 
the dynamics of cell polarisation. In these models, molecular markers, such 
as proteins, diffuse in the cytoplasm and are actively transported along the 
cytoskeleton. The resulting motion is a biased diffusion regulated by the 
markers themselves. Using numerical simulations and mathematical heuris- 
tics, we observe that the coupling on the velocity field achieves an inhomo- 
geneous distribution of molecular markers without any external asymmetric 
field. Such an inhomogeneous distribution is only due to interaction between 
molecular markers. 

Throughout this paper, the density of molecular markers (resp. advec- 
tion field) is denoted by p(t,x) (resp. u(t,x)). The advection is obtained 
though a coupling with the membrane concentration of markers. The cell is 
figured by the domain £1 C M n with n = 1 , 2 and a part of the boundary of 
the domain will be the active membrane denoted by T. The time evolution 
of the molecular markers satisfies the following advection-diffusion equation, 
see [H] and [3]: 



There is no creation nor degradation of molecular markers in the cell, so the 
quantity of molecular markers remains constant in time: 



d t p(t,x) = DAp(t,x)-xV. (p(t,x)u(t,x)), t>0, xe!l, 
p(0,x) = /3 (x). 



(1) 
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This condition is ensured by a zero flux boundary condition on the boundary. 
A first simplified step is to assume that the cell is essentially bidimensional 
and to neglect curvature effects. The membrane boundary is then a ID line 
along the y-axis and the cytoplasm is parametrized by x = (x, y) G R+ x R. 

The plan of this work is the following. First, we recall the main mathe- 
matical results of the simplified model in ID for 0, = (0, oo) and T = {x = 0}, 
see [1], [3] for more details. Then we study a more realistic model, that in- 
cludes dynamical exchange of markers on the boundary for a general f2. 
This model was introduced in [8] and studied in [3] in the one dimensional 
case. Here, we will perform a first numerical analysis of this model in the 
two dimensional case, for periodic (in one direction) and bounded (in the 
other direction) domain. Finally, we provide a methodology for parameter 
estimation by using mathematical heuristics and biological literature. 

1.1 One dimensional case 

In this section, we study the one dimensional case on the half line for f2 = 
(0, oo). The membrane is then the point T = {x = 0}. For the first model, 
the advection field towards the membrane is equal to the density of molecular 
markers on the boundary p(t, 0). Then we improve this model by considering 
that only the trapped molecular markers on the membrane contribute to the 
advection field. 

1.1.1 Simplified model set on the half line 

In [3] a first mathematical studies has been done on this model. We define 
an advection field u(t, x) for ([!]) 

u(t,x) = -p(i,0), 

in such a case reads as (with D = 1 and x = 1) : 

d t p(t,x) = d xx p(t,x) + p(t,0) d x p(t,x), t>0, x > 0, (3) 

with the following zero flux condition on the boundary {x = 0}, that ensures 
the mass conversation Q, 

d xP (t, 0) + p{t, Of = 0. (4) 

In [3] , it has been proved that solutions of ^ blow-up in finite time if their 
masses are above a certain critical mass, M > 1, and exist globally in time 
if M < 1. Let us first recall the definition of weak solutions of ([3]). 

Definition 1.1 We say that p(t,x) is a weak solution of Q on (0,T) if it 
satisfies: 

P €L co (0,T-L 1 + (R + )), d x peL l ((0,T)xR + ), 
and p(t,x) is a solution of (p| in the sense of distributions in V(H + ). 
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Let us now recall the main results for weak solutions of (13 



Theorem 1.2 (Global existence: M < 1) Assume that the initial data 
po satisfies both po G L l {(\ + x) dx) and f x>Q po(x)(log po(x)) + dx < +oo. 
Assume in addition that M < 1, t/ien there exists a global weak solution of 
equation 

Theorem 1.3 (Blow-up: M > 1) Assume M > 1. ^4ny weak solution of 
equation ^ with non-increasing initial data po blows-up in finite time. 

Remark 1 It would tempted to interpret blow-up of solutions of the one 
dimensional model as cell polarisation. But it is to be noticed that con- 
centration of markers on the boundary doesn't mean polarisation. Indeed, 
consider a radially symmetric 2D cell case. Equation then reduces to the 
one dimensional one. Above a threshold on the total mass, the convection 
wins and markers concentrate on the boundary. In some situations, these 
markers may be homogeneously distributed on the boundary and in such a 
case there is no symmetry breaking. 

1.1.2 The model with dynamical exchange of markers at the 
boundary 

Such a direct activation of transport on the boundary seems to be unrealistic. 
Indeed possible occurrence of blow-up in finite time suggests this claim. We 
improve the previous model by distinguishing between cytoplasmic content 
p(t, x) and the concentration of trapped molecules on the boundary, that will 
be denoted by p(t). The dynamical exchange of markers at the boundary 
is done with an attachment rate k on and a detachment rate k ff, hence the 
time evolution of p(t) is 

^{t) = kon p{t, 0) - k ff p{t). (5) 
The advection field u(i, x) in ([!]) is now defined by 

u(t,x) = -p(t), 
hence Q (with D = 1 and x = 1) reads as: 

d t p(t,x) = d xx p(t,x) + p(t)d x p(t,x), t>0, x>0, (6) 
with a modified boundary condition 

d x p(t,0) + p(t,0)p(t) = ^-u(t). (7) 
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This ensures the following mass conservation shared among p(t, x) and p(t): 



M = I po(x)dx + /xo = / p(t,x)dx + p(t). (8) 
Jr+ ' Jr + 

With equation the self-activation of transport by p(t, 0) is then delayed 
in time. Since the transport speed is bounded p(t) < M, the solution of the 
model with dynamical exchange on the boundary exists globally in time. 
More precisely it is possible, see [3], to prove that it converges towards a 
non trivial stationary state. 



1.2 Two dimensional case : the model with dynamical ex- 
change of markers at the boundary 

Let f2 C M. 2 be the cytoplasm domain, as in the one dimensional case ([5]) we 
consider dynamical exchange of markers at the boundary, so for x G T we 
have the evolution in time of p(t, x) 

d t p(t, x) = kon p(t, x) - k ff fj,(t, x) . (9) 

with a modified boundary condition for p(t, x) at point x G T 

(DVp(t,x) - X p(t,x)u(t,x)).n x = -%<t,x). (10) 

where n x is the outward normal to T. This ensures the following mass 
conservation sharing by /f(i,x) and p(t,x): 

M= po( x )^ x + / W)( x )^ x = / p(i,x)<ix+ / /i(t,x)cbc. (11) 
Jn Jr Jn Jr 

We consider the following advection field deriving from a harmonic potential 
modeling the transport by actin filaments (cytoskeleton): 



u(i, x) = Vc(t,x), where 



-Ac(i,x) = 0, 



if x G Q, 



Vc(t,x).n x = 5(x)/i(t,x), if x G T. 
This advection field orientation is due to the actin networks. 



(12) 




Actin filaments are attached on the membrane and randomly distributed, 
there orientations are mixed up. We also add the external pheromone con- 
centration at x G r which acts by the mating-pheromone MAPK cascade 
on the actin transport. 

In dimension 2, we have global existence for the model without exchange 
on the boundary (replacing equation (12) by Vc(i,x).n x = S(x)p(t,x) if x G 
r) with f2 = (0, +00) x E and V = {0} x E. For clarity, we recall this result, 
sec [3j for more details. 

Theorem 1.4 (Global existence in dimension 2) Assume that the ad- 
vection field satisfies the two following conditions: V • u > and u(t, 0, y) ■ 
e e = p(t,0,y). Assume that the initial data po satisfies both po G -^ 1 ((1 + 
|x| 2 )dx) and \\po\\l 2 smaller than some constant c. Then there exists a 
global weak solution to equations Q-Q. 



In the two dimensional case, for the model with exchange on the bound- 
ary, blow-up or global existence have not been proved yet. In this work, 
we make a first step in this direction by using a mathematical heuristic and 
numerical simulations. 



1.3 Heuristics 

The mathematical analysis performed in [3] has demonstrated that a class 
of models exhibit pattern formation (either blow-up or convergence towards 
a non homogeneous steady state) under some conditions. However the main 
question still remains unanswered: do these models describe cell polarisation 
or not? Thus in order to provide a first answer to this question, we will 
perform numerical simulations. Our aim is to see if, under some conditions, 
the model leads to a concentration of markers, not only on the boundary, 
but on a small region of the boundary. In such Si CELS6 polarisation occurs. 
In order to obtain more information on the critical value distinguishing 
the polarised case and the stable case, in the two dimensional case we will 
use a mathematical heuristics that we describe now. Let x = (x, y) be in 
ft = R + x R and let V = {0} x R be the boundary, we have that 



u(t, x) = Vc(t, x), where 



-Ac(t,x) = 0, ifxG 
-d x c(t,0,y) = S(y)/j,{t,y), 'dye'. 



hence, see [7] e.g., it is well known that 

c{x,y) = --[ log( \/ (y - y') 2 + x 2 )(Sp)(y')dy' . 

n Jy'eR 

The tangential component at the boundary is then given by 
u(t,0,y)-e y = -H(Sp){y), yeR, 



(13) 
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where H denotes the one-dimensional Hilbert transform that we recall now, 
see [5] e.g., with respect to the y variable: 

Wfj)(y) = -p.v. / -^f(y')dy'. 
Jrv-v' 

Integrating the main equation ([I]) with respect to x with zero flux condition 
on r = {x = 0}, we obtain: 



d t / p(t,x,y)dx = Ddyy I / p(t, x, y) dx J ~x d y 
Jx>0 \Jx>0 J 



p(t,x,y)(u(t,x,y) ■ e y )dx 

x>0 



In the super-critical case, numerical simulations, see [IT], suggest that the 
density p(t, x) concentrates on the boundary {x = 0}. Assuming p(t, x, y) = 
is(t,y)5(x = 0), we can formally write the dynamics of v{t,y) as follows: 

d t u{t,y) =Dd yy u(t,y) + xdy(u(t,y)H(Sp J )(y)) . 

Assuming S constant on R and p(t,y) = Yjjv(t,y) for y G R, it reads as 

d t u(t, y) = D d yy v{t, y) + d y (u(t, y)U{v){y)) . 

k off 

Hilbert transform has a critical singularity to offset the diffusion on this 
equation [5]. We have a blow-up if J R v(t, y) dy = M is above 2 g^° JJ . This 
is a first step to observe a critical mass phenomenon and this may lead to 
blow-up if the mass is large enough. In this way, we define an order of 
magnitude for some parameters. 

It is to be noticed that this latter criterion is valid for an infinite domain, 
namely y £ M. In the case of a cell, the domain will be finite and the 
existence of such a dichotomy has not been proved yet. In order to see if 
such a dichotomy holds true we will perform numerical simulations. This is 
the object of the following section. 



2 Numerical analysis 

We first give a discretization of the convection-diffusion model set on a ID 
periodic domain. This first step allows us introducing the discretization of 
this model on a 2D domain which is periodic in one direction and bounded 
on the other direction. 



2.1 One dimensional case 

Let u(t, x) be a given function. We consider the following advection-diffusion 
equation on the periodic domain £1 = IR/Z 

d t p{t,x) = d x (d x p(t,x) - u(t,x) p(t,x)), t>0, xeQ. (14) 



7 



Let t n = ndt be the time discretization and {xj = jdx,j £ {1, ...,N X }} be 
the space discretization of the periodic interval IR/Z. Since the equations of 
the model are written in a conservative form, the natural framework to be 
used for the spatial discretization is the finite volume framework. We hence 
introduce the control volume defined for j £ {!,•••, N x } 



V; 



(15) 



Let (resp. ) De the approximated value of the exact solution p(t n , Xj) 



(resp. u(t n , x ■ , i )), the classical upwind scheme for (14) reads as 

J ' 2 ■ ' 



J 2 



dt dx ' ^ {!,...,*,}, 

where the numerical flux T- , i and .F-_i are defined by 



(16) 



J 2 



n+1 



„n+l 



A up {u 



„n+l 



n+1 
^-1 



with the advection numerical flux is given by 



^(«?_l,Pi-i,#). 



A up (w,x_,x+) 




si ti > 0, 
si u < 0. 



(17) 



The periodic flux condition on boundary reads as T\ = J r Nx+ i and we set 
the value it" = u 1 ^ 1 . The diffusion part is treated implicitly and it is then 

unconditionally stable, while the advection term is treated explicitly. The 
CFL condition of the scheme is 



< 



dx 
~dt' 



We define the column vector p n = (p™ pV^ ... Pn ) ■ As usual, see e.g. 
p], the discrete heat matrix A S Mn x (M) with periodic flux condition on 
the boundary is defined as 



/n _i_ dx 2 
• ' dt 



A 



-1 2 + 



1 

dx 2 
dt 



\ -i 



• 2 + 



dx 2 
dt 



"I \ 



9 _l_ ^ 2 / 

dt f 



(18) 
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Periodic flux condition adds the top right term and the bottom left term. 
Next, in order to use A up defined by equation (17), we define {u) + = 
max(u, 0) and (u)~ = min(u, 0). The discrete advection matrix B G Mn x (M) 
with periodic flux condition on the boundary is then defined as in [1] 



D 



dx' 
lit 



+dx 



u'l 



u n 3 



V 



In -dx 



N x - 



J 2 



u 



'3+1 



u n 1 

J 2 



u 



N x - 



u 



'3+1 



U 



N. r , 



II 



I 



(19) 



We use a standard numerical method to invert the symmetric positive defi- 
nite matrix A. Finally, at each time step we resolve 



A~ l Bp n . 



2.2 Two dimensional case 

We perform numerical simulations on the model with dynamical exchange 
of markers at the boundary. In this work, we assume that the cell occupies 
a disk of radius r > 0. Furthermore for simplicity, we consider a bounded- 
periodic domain = [0, r] x M/27rrZ with T = {r} x ]R/27rrZ. This simplifies 
our numerical approach by using finite difference schemes on Cartesian grid. 
We start with the numerical study of the equation on p by assuming that 
the advection field u(t, x) = Vc(i,x) is known. Then we perform the dis- 
cretization of c. 

In this section, for simplicity we fix all the parameters values to 1 except 
M. Let us first recall the model with dynamical exchange of markers at the 
boundary on = [0, r] x M/27rrZ: 

8 tP = V. (Vp - p Vc) , in (0, r) x M/2vrrZ, (20) 
d x p - pd x c = -d t p, on {r} x M/27rrZ, (21) 
d x p - pd x c = 0, on {0} x R/2vrrZ. (22) 
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Dynamical exchange markers on active boundary {r} x IR/27rrZ is given by 

d t p = p- p, on {r} x M/2vrrZ. (23) 

Laplace equation on c with non appropriate Neumann conditions on a bounded 
domain is ill-posed, see PQ e.g. In order to handle this problem, we add a 
degradation term: 



Ac + a c = 0, in (0, r) x R/2irrZ, 
—d x c = p, on {r} x R/2irrZ, 
-0 x c = 0, on {0} x M/2vrrZ. 



(24) 
(25) 
(26) 



We take random initial conditions c, po an d po satisfying the following mass 
conservation 



Po + / P-o 



M. 



(27) 



Let i™ = ndt be the time discretization and {xj = jdx,j £ {1, A^}} 
(resp. {yk = kdy,k G {1, N y }}) be the space discretization of the 
bounded interval [0,r) (resp. periodic interval M/27rrZ). We introduce the 



control volume Wi 



W, 



C 



(?,*) 



(^_i,^+0 x (y fc _i,y fe+ i). 



Let (resp. p 7 ^) be the approximated value of the exact solution 



p(t n ,Xj,yk) of equations (20), (21), (22) and (27) (resp. p(t n ,yk) of equations 



p3[ ) and (27)). Let C (j t k) be the approximated value of the exact solution 
c{xj,yk) of equations (24), (25) and (26). 



2.2.1 Equation on p 

We can resolve at each time step for k G {1, N y } 

p^=pl + dt(pl-pl). 

2.2.2 Equation on c 

For simplicity, we call J- the numerical flux as in the ID case, we can write 
the following scheme: for (J, k) G {1, A^} x {1, N y } 



dx 



+ 



with numerical flux defined by 



dy 



c (j+l,k) ~ C(j,k) 



ac U,k) = o. 



(j+hk) dx 

c (i,fc) ~ c (i-i,fc) 
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The zero flux boundary condition (26) impose that Tn k \ 
boundary condition (25) F (N , 1 



and the 



for k G {1, N y }. Similarly, the 
periodic conditions impose T,- N ,i-v = J-j • u for j G {1, ...,N X }. We define 



the column vector C by C(k + (j — l)N y ) = Cuja with (j, k) G {1, A^.} x 
{1, ...,N y }. As previously the rigidity matrix A^d^ol is defined by 



( A a + Id -Id 
-Id A a + 2Id 

-Id 



A 



2D,a 



-Id 

A a + 2Id 



\ 



A a + 2Id -Id 

-Id A a + 2Id -Id 

-Id A a + IdJ 



where discrete Poisson matrix A a G M^ y (M) in ID with periodic flux con- 
dition on boundary is defined by 



A, 



(2 + adx 2 
-1 



V 



-1 
2 + adx 2 



\ 



2 + adx 2 -1 

-1 2 + adx 2 J 



The flux boundary condition {r} x M/27rrZ imposes this right hand side 
column vector of length N x N y : 



R, 



-dx((^) k 



0). 



We use a standard numerical method to invert the symmetric positive defi- 
nite matrix Aid a and then resolve at each time step 



C 



2.2.3 Equation on p 

For simplicity, we call J- the numerical flux as in the previous cases, we can 
write the upwind scheme as follows: 

dt dx dy 
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where = 
flux is denned by 



(j+l,fc) U,k) „.n 



dx 



dx 



and the numerical 



n+l _ n+l 

n+ l _ n+ l 

da; 



U+h,ky 0'.*)' (i+MOJ ' 



The zero flux boundary conditions (22) impose Jvi fc % 



0, while the flux 



boundary conditions (j 2 1 1) impose T/ Nx+ i k \ = ~^~dT~ M ~ ^ or ^ e {!> ••> ^i/}- 

T, ,ix for j G {1,...,JV X }. 



Similarly, the periodic conditions impose J 7 , 



We define the column vector V n by 7 :>n (/c + (j — l)N y ) = PJ\ ^ with (j, k) G 
{1,...,^} x {1, iVy}. In what follows for simplicity we consider that 
dx = dy. We define the rigidity matrix A 2 d 6 M NxNy (R) with A G M Ny (R) 
defined by (pi): 



-Id A + 2/d 



2D 



A + 2Id -Id 
-Id A + Id J 



We define the following diagonal matrices for j G {1, N x }, UT^x G 
Mjv b (R) and 17" , G M Ny (R): 

2 



V 

(■ 



U7 ! 

J+3 



V 



(i+|,fc-l) J 



(*+§,*)' 



(11 



•7 



With 5 G Mtv (M) defined by equation (19), the discrete advection matrix 



B2D £ Mjy x N y (M) with zero flux boundary condition in the x-axis direction 
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and periodic flux boundary condition in the y-axis direction is defined by 



/Ut U 



B 



2D 



B 



B 



dx 



u 
u 



(VI 



+ dx 



U 



Ul 



u: 



\ 



K x u 



N x - 



The flux boundary condition {r} x M/27rrZ imposes this right hand side 
column vector of length N x N y : 



Rc 



-dx((^Z^) fc o ... o)- 



We use a standard numerical method to inverse the symmetric positive def- 
inite matrix A^d and then resolve at each time step 

V n+1 = A 2 '(B 2D V n + R p ). 



2.3 Graphics 

With the previous numerical analysis, we implement all numerical simula- 
tions using MATLAB. We test different values of M: 



3 Conclusion 

In this work we have provided a first answer to the following question: do 
the nonlinear convection-diffusion models given in [H] and [3] describe cell 
polarisation or not? To do so we have used both a mathematical heuristic 
and numerical simulations. Numerical simulations were necessary because 
the heuristic is only valid for an infinite geometry while the cell is obviously 
finite. The numerical simulations ensure that solutions develop symmetry 
breaking over a critical value M* given us a first justification of the mathe- 
matical heuristic. In a further work, we will estimate an approximate value 
of this critical mass. 
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A) B) 




Figure 1: Numerical Simulations on 17 = [0, 1] x M/27rZ with T = {r} x 
M/27rrZ and all parameters equal to 1. A) For M = 20 greater enough, 
symmetry breaking appears. Molecular markers are concentrated on one 
point of the membrane in finite time. B) For M = 0.01 small, steady state 
is homogeneous in the y-axis. 
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